﻿using IntelligentEscalator.Models;
using IntelligentEscalator.Tools;
using System;
using System.Collections.Generic;
using System.Numerics;
using System.Reflection.Metadata;
using System.Runtime.InteropServices;
using System.Text;
using System.Threading;

namespace IntelligentEscalator.DLL
{
    public class FFTWWrapper
    {
        private static readonly object _lock = new object();    
        private static Mutex mutex = new Mutex();
        public FFTWWrapper()
        {

        }

        // 计算FFT频谱
        public ResultData CalculateFftSpectrum(double[] filedata, int Fs, int N)
        {
            lock (_lock)
            {
                //mutex.WaitOne();
                int length = 2 * N;
                ResultData resultData = new ResultData();
                QCPGraphData[] fftmDataArray = new QCPGraphData[N / 2];

                IntPtr mfin = fftw.malloc(length * 8);
                IntPtr mfout = fftw.malloc(length * 8);
                // 分配空间
                double[] din = new double[length];
                double[] dout = new double[length];
                for (int i = 0; i < N; i++)
                {
                    din[2 * i] = filedata[i];
                    din[2 * i + 1] = 0;
                }
                // 将输入信号填充到复数数组中
                fftw.SetData(mfin, din, length);
                //Marshal.Copy(din, 0, mfin, length);
                // 创建 FFT 计划
                IntPtr plan = fftw.dft_1d(N, mfin, mfout, fftw_direction.Forward, fftw_flags.Estimate);


                // 执行 FFT
                fftw.execute(plan);



                double[] spectrummagnitude = new double[N / 2];
                //double[] complex = mfout.GetData_Double();
                double[] complex = new double[length];
                Marshal.Copy(mfout, complex, 0, length);
                complex = fftw.GetData(mfout, length);

                //Complex[] complex = mfout.GetData_Complex();
                // FFT频谱计算
                for (int i = 0; i < N / 2; i++)
                {
                    double magnitude = Math.Sqrt(complex[2 * i] * complex[2 * i] + complex[2 * i + 1] * complex[2 * i + 1]);
                    //double magnitude = Math.Sqrt(complex[i].Real * complex[i].Real + complex[i].Imaginary * complex[i].Imaginary);
                    spectrummagnitude[i] = magnitude / N * 2;

                    var fftNewPoint = new QCPGraphData
                    {
                        Key = (double)i * Fs / N,
                        Value = magnitude / N * 2,
                    };
                    fftmDataArray[i] = fftNewPoint;
                }
                // 销毁计划 释放内存
                fftw.destroy_plan(plan);
                fftw.free(mfin);
                fftw.free(mfout);

                //mutex.ReleaseMutex();
                resultData.qCPGraphDatas = fftmDataArray;
                resultData.data = spectrummagnitude;
                return resultData;
            }
        }

        // 倒谱计算
        public QCPGraphData[] calculate_Cepstrum(double[] Spectrummagnitude, int Fs, int N)
        {
            double[] logSpectrum = new double[N / 2]; // 长度是 N / 2
            QCPGraphData[] cepstrumDataArray = new QCPGraphData[N / 2];

            IntPtr ifft_in = fftw.malloc(N / 2 * 16);
            IntPtr ifft_out = fftw.malloc(N / 2 * 16);

            // 对数谱
            for (int i = 0; i < N / 2; i++)
            {
                logSpectrum[i] = Math.Log(Spectrummagnitude[i] + 1e-10); // 加小常数以避免对数为负无穷
            }

            logSpectrum = RemoveMean(logSpectrum);   // 去均值 直流分量

            // 执行ifft逆变换 
            // 反变换到频域
            double[] i2 = new double[N];
            for (int i = 0; i < N / 2; i++)
            {
                i2[2 * i] = logSpectrum[i];
                i2[2 * i + 1] = 0;
            }
            // 将输入信号填充到复数数组中
            fftw.SetData(ifft_in,i2,N);
            // 执行ifft逆变换
            IntPtr plan = fftw.dft_1d(N/2,ifft_in,ifft_out, fftw_direction.Backward, fftw_flags.Estimate);
            // 执行逆FFT
            fftw.execute(plan);

            fftw.destroy_plan(plan);
            double[] complex = fftw.GetData(ifft_out, N); 
            // 计算倒谱 ;
            for (int i = 0; i < N / 2; i++)
            {
                double magnitude = Math.Sqrt(complex[2 * i] * complex[2 * i] + complex[2 * i + 1] * complex[2 * i + 1]);
                //Spectrummagnitude[i] = magnitude / N * 2;
                var fftnewPoint = new QCPGraphData
                {
                    Key = (double)i * Fs / N,
                    Value = magnitude / N * 2,
                };

                cepstrumDataArray[i] = fftnewPoint;

            }

            //释放fft内存
            fftw.free(ifft_in);
            fftw.free(ifft_out);
            return cepstrumDataArray;

        }
        // 包络谱
        public QCPGraphData[] calculate_envelope_spectrum(double[] filedata, int Fs, int N)
        {
            // 分配空间
            IntPtr  env_in = fftw.malloc(N * 16);
            IntPtr env_out = fftw.malloc(N * 16);

            double[] spectrum = new double[N];
            double[] envelope = new double[N];
            QCPGraphData fftnewPoint;
            QCPGraphData[] spectrumDataArray = new QCPGraphData[N / 2];

            double[] i2 = new double[2 * N];
            for (int i = 0; i < N; i++)
            {
                i2[2 * i] = filedata[i];
                i2[2 * i + 1] = 0;
            }
            fftw.SetData(env_in, i2,2*N);
            IntPtr plan = fftw.dft_1d(N, env_in, env_out, fftw_direction.Forward, fftw_flags.Estimate);
            fftw.execute(plan);
            fftw.destroy_plan(plan);

            // 实现 Hilbert 变换
            double[] i3_out = fftw.GetData(env_out,2*N);
            double[] o3 = new double[2 * N];
            for (int i = 0; i < N; i++)
            {
                if (i > 0 && i < N / 2)
                {
                    o3[2 * i] = i3_out[2 * i] * 2;
                    o3[2 * i + 1] = i3_out[2 * i + 1] * 2;
                }
                else if (i == 0 || i >= N / 2)
                {
                    o3[2 * i] = 0;
                    o3[2 * i + 1] = 0;

                }
            }
            fftw.SetData(env_out,o3,2*N);
            // 执行逆 FFT
            plan = fftw.dft_1d(N, env_in, env_out, fftw_direction.Backward, fftw_flags.Estimate);
            fftw.execute(plan);

            fftw.destroy_plan(plan);
            double[] complex =fftw.GetData(env_in,2*N);
            // 计算包络
            for (int i = 0; i < N; i++)
            {
                double magnitude = Math.Sqrt(complex[2 * i] * complex[2 * i] + complex[2 * i + 1] * complex[2 * i + 1]); // 计算幅度
                envelope[i] = magnitude / N * 2;
            }

            envelope = RemoveMean(envelope);   // 去均值 直流分量
            // 清理

            //-------------- 计算包络谱
            double[] eni = new double[2 * N];
            for (int i = 0; i < N; i++)
            {
                eni[2 * i] = envelope[i];
                eni[2 * i + 1] = 0;
            }
            fftw.SetData(env_in,eni,2*N);
            // 执行 FFT
            plan = fftw.dft_1d(N, env_in, env_out, fftw_direction.Forward, fftw_flags.Estimate);
            fftw.execute(plan);


            complex = fftw.GetData(env_out, 2 * N);
            // 计算包络谱
            for (int i = 0; i < N; i++)
            {
                double magnitude = Math.Sqrt(complex[2 * i] * complex[2 * i] + complex[2 * i + 1] * complex[2 * i + 1]); // 计算幅度谱
                spectrum[i] = magnitude / N * 2;
            }

            // 计算包络谱
            for (int i = 0; i < (N / 2); i++)
            {
                fftnewPoint = new QCPGraphData();
                fftnewPoint.Key = (double)i * Fs / N;
                fftnewPoint.Value = spectrum[i];
                spectrumDataArray[i] = fftnewPoint;
            }


            // 清理
            fftw.destroy_plan(plan);
            fftw.free(env_in);
            fftw.free(env_out);
            return spectrumDataArray;
        }
        // 功率谱 
        public QCPGraphData[] calculate_power_spectrum(QCPGraphData[] fftmDataArray, int Fs, int N)
        {
            QCPGraphData[] powerSpectrumDataArray = new QCPGraphData[N / 2];
            QCPGraphData fftnewPoint;

            for (int i = 0; i < N / 2; ++i)
            {
                fftnewPoint = new QCPGraphData();
                fftnewPoint.Key = (double)i * Fs / N;
                fftnewPoint.Value = fftmDataArray[i].Value * fftmDataArray[i].Value;   // 功率谱
                powerSpectrumDataArray[i] = fftnewPoint;
            }
            return powerSpectrumDataArray;
        }
        // 特征谱
        public FeatureModel calculate_feature_spectrum(double[] filedata)
        {
            double max_value = 0;
            double min_value = 0;
            double abs_max_value = 0;

            // 最大值 最小值 峰值计算
            for (int i = 0; i < filedata.Length; i++)
            {
                if (filedata[i] > max_value)
                {
                    max_value = filedata[i];
                }
                if (filedata[i] < min_value)
                {
                    min_value = filedata[i];
                }

                if (Math.Abs(filedata[i]) > abs_max_value)
                {
                    abs_max_value = Math.Abs(filedata[i]);
                }
            }
            FeatureModel featureModel = new FeatureModel();
            featureModel.acc_peakToPeak = (float)(max_value - min_value);                                                                        // 峰峰值 OK
            featureModel.acc_peak = (float)(abs_max_value);                                                                                // 峰值 OK
            featureModel.acc_meanAbsolute = (float)calculateMeanAbsolute(filedata, filedata.Length);      // 绝对均值 OK
            featureModel.acc_rmsValue = (float)calculateRMS(filedata, filedata.Length);               // 有效值（均方根值） OK
            featureModel.acc_kurtosisValue = (float)calculateKurtosis(filedata);                                                   // 峭度  OK
            featureModel.acc_Wave = (float)(featureModel.acc_rmsValue / featureModel.acc_meanAbsolute);                      // 波形指标  OK
            featureModel.acc_peakFactor = (float)(featureModel.acc_peak / featureModel.acc_rmsValue);                              // 峰值指标  OK
            featureModel.acc_PulseFactor = (float)(featureModel.acc_peak / featureModel.acc_meanAbsolute);                          // 脉冲指标  OK
            featureModel.acc_CrestFactor = (float)(featureModel.acc_peak / calculateaqrtRMS(filedata, filedata.Length));  // 裕度指标
            featureModel.acc_kurtosisFactor = (float)(featureModel.acc_kurtosisValue / Math.Pow(featureModel.acc_rmsValue, 4));             // 峭度指标  OK
            featureModel.acc_envelope = 0;     // 包络
            featureModel.acc_rmsvelocity = 0;
            return featureModel;
        }

        #region
        // 计算绝对均值
        public double calculateMeanAbsolute(double[] data, int size)
        {
            double sum = 0.0;
            for (int i = 0; i < size; i++)
            {
                sum += Math.Abs(data[i]);
            }
            double meanAbsolute = sum / size;
            return meanAbsolute;
        }
        // 计算均方根值
        public double calculateRMS(double[] data, int size)
        {
            double sum = 0.0;
            for (int i = 0; i < size; i++)
            {
                sum += data[i] * data[i];
            }
            double meanSquare = sum / size;
            double rms = Math.Sqrt(meanSquare);
            return rms;
        }
        // 去除平均值
        public double[] RemoveMean(double[] data)
        {
            double sum = 0.0;
            double mean = 0.0;
            double[] dataFilter = new double[data.Length];
            for (int i = 0; i < data.Length; i++)
            {

                sum += data[i];
            }

            mean = sum / data.Length;

            for (int i = 0; i < data.Length; i++)
            {

                dataFilter[i] = data[i] - mean;
            }
            return dataFilter;
        }
        // 计算峭度
        public double calculateKurtosis(double[] data)
        {
            double sum = 0.0;
            for (int i = 0; i < data.Length; ++i)
            {
                sum += Math.Pow(data[i], 4);
            }
            return (sum / data.Length);
        }
        public double calculateaqrtRMS(double[] data, int size)
        {
            double sum = 0.0;
            for (int i = 0; i < size; i++)
            {
                sum += Math.Sqrt(Math.Abs(data[i]));
            }
            double meanSquare = sum / size;
            double rms = meanSquare * meanSquare;
            return rms;
        }

        // 计算平均值
        public double calculateMean(double[] data)
        {
            double sum = 0.0;
            for (int i = 0; i < data.Length; ++i)
            {
                sum += data[i];
            }
            return sum / data.Length;
        }
        #endregion

        public ChannelData CalculateChannelData(int channel,double[] filedata, int Fs, int N) {
            filedata = RemoveMean(filedata);
            try {
                ChannelData channelData = new ChannelData();
                ResultData result = CalculateFftSpectrum(filedata, Fs, N);
                //QCPGraphData[] ifftw = calculate_Cepstrum(result.data, Fs, N);
                //QCPGraphData[] envData = calculate_envelope_spectrum(filedata, Fs, N);
                //QCPGraphData[] power = calculate_power_spectrum(result.qCPGraphDatas, Fs, N);
                //FeatureModel featureModel = calculate_feature_spectrum(filedata);
                channelData.ChannelIndex = channel;
                channelData.Datas = filedata; // 原始数据
                channelData.FFTDatas = result.qCPGraphDatas; // FFT数据
                //channelData.IFFTDatas = ifftw; // 倒谱数据
                //channelData.ENVDatas = envData; // 包络谱数据
                //channelData.PownerDatas = power; // 功率谱数据
                //channelData.FeatureModels = featureModel; // 特征值
                return channelData;
            }
            catch (Exception e) {
                Console.WriteLine("数据转换:"+e.ToString());
                LogHelper.Error("数据转换:"+e.ToString());
                return null;
            }
        
        }
    }


    public class QCPGraphData
    {
        public double Key { get; set; }
        public double Value { get; set; }
    }
    public class ResultData
    {
        public double[] data { get; set; }
        public QCPGraphData[] qCPGraphDatas { get; set; }
    }
}
